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ABSTRACT 

We study the formation and long-term evolution of primordial protostellar disks harbored by first 
stars using numerical hydrodynamics simulations in the thin-disk limit. The initial conditions are 
specified by pre-stellar cores with distinct mass, angular momentum, and temperature. This allows 
us to probe several tens of thousand years of the disk's initial evolution, during which we observe 
multiple episodes of fragmentation leading to the formation of gravitationally bound gaseous clumps 
within spiral arms. These fragments arc torqucd inward due to gravitational interaction with the spiral 
arms on timescales of 10 3 -10 4 yr and accreted onto the growing protostar, giving rise to accretion and 
luminosity bursts. The burst phenomenon is fueled by continuing accretion of material falling onto the 
disk from the collapsing parent core, which replenishes the mass lost by the disk due to accretion, and 
triggers repetitive episodes of disk fragmentation. We show that the burst phenomenon is expected 
to occur for a wide spectrum of initial conditions in primordial pre-stellar cores and speculate on how 
the intense luminosities (~10 L©) produced by this mechanism may have important consequences 
for the disk evolution and subsequent growth of the protostar. 

Subject headings: cosmology: theory, stars: formation, accretion disks, hydrodynamics 



1. INTRODUCTION 

Cosmological-scale simulations of collapsing primor- 
dial clouds in the early universe have been used to sug- 
gest that the first luminous objects in the universe were 
stars with masses of M > 100 M^ that formed in rel- 
ative isolation (e g ., lAbell l2000t iBromm k Loebl 12004 
lYoshida et al.|[2008 ). However, the above numerical sim- 
ulations were effectively based on a scenario of monolithic 
quasi-spherical collapse of cloud cores. Angular mo- 
mentum, disk formation, and fragmentation have added 
greater depth to this picture, replacing it with one in 
which the collapsing primordial cores of the early uni- 
verse produce rich structure in the inner regions where 
disks emerge. 

Although protostellar disks are a ubiquitous outcome 
of the present-day star formation process, the importance 
of these structures to the evolution o f the first stars has 
beg un to be understood on ly recently. ISaigo et all (j2004l ) 
and iMachida et al.l (|2008| ) have demonstrated that disk- 
like structures are expected to form from primordial cores 
with a wide range of initial rotation rates, a nd can frag- 
ment t o yield binary pairs of first stars. IClark et al.l 
were able to follow the collapse of primordial 
clouds to protostellar densities and found vigourous grav- 
itational fragmentation in primordial disks leading to the 
formation of tightly bound multiple stellar systems. Us- 
ing smoothed particle hydrodynamics simu lations com- 
bined with the sink particle technique, ISmith et al.l 
(|2011l 120121 ) studied the importance of accretion lumi- 
nosity on disk fragmentation and the effect of protostel- 
lar ac cretion on the str ucture and evolution of proto- 
stars. iGreif et "all (|2012f ) most recently performed high- 
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resolution, three-dimensional hydrodynamic simulations 
using a Lagrangian moving mesh and found rapid migra- 
tion and merging of secondary fragments with the pri- 
mary protostar. 

The general picture to emerge from these simula- 
tions is that protostellar accretion may be a highly vari- 
able process fuelled by disk gravitational instability and 
fragmentat ion, in a similar manner to present-day star 
formation (jVorobvov k Basull2006l 120101 : IMachida et "all 
120111 ). However, a major problem of high-resolution 
three-dimensional studies has been the difficulty of fol- 
lowing the evolution of primordial protostellar disks for 
longer than a few thousand years. In this study we 
present gas hydrodynamics simulations of collapsing pri- 
mordial cores into a disk formation phase using the thin- 
disk approximat ion. We adopt a bar otropic equation of 
state derived bv lOmukai et al.1 ()2005l ) for the primordial 
chemical composition of gas. The benefit of this type 
of formulation is that we are able to study the fragmen- 
tation, evolution and subsequent accretion of protostel- 
lar mass clumps formed within the disk, while simulta- 
neously maintaining resolution of the extended remnant 
parent cloud core over many orbital periods and model 
realizations. 

The structure of this article is as follows. We briefly 
describe the numerical simulations in Section 2, includ- 
ing the modifications that have been made pertaining to 
the primordial star-forming environment. In Section 3 
we present characteristics of the temporal evolution of 
our reference model, starting from the prestellar phase, 
and ending when the mass of the protostar has reached 
45 Mq. Section 3 also includes an analysis of the pro- 
tostellar accretion luminosity expected to arise from the 
burst mode of accretion that is unique to our model of 
primordial star formation. Finally, in Section 4 we ex- 
tend our discussion and draw conclusions from these re- 
sults. 

2. MODEL DESCRIPTION 
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Our model and method of so l ution update the model 
presented in IVorobvov fc Basul (|2005al l2006| ) , with ap- 
propriate modifications for star and disk formation in the 
early universe. We follow the evolution of gravitationally 
unstable primordial cores from the isolated pre-stellar 
stage into the protostar and disk formation stage and 
terminate our simulations once about 30% of the initial 
mass reservoir has been accreted onto the protostar plus 
disk system. Once the disk is formed, it occupies the in- 
nermost region of our numerical grid, while the infalling 
envelope — the remnant of the parent core — occupies the 
outer extent. The dynamics of both the disk and enve- 
lope are followed self-consistently on one global numerical 
grid space, which ensures correct mass infall rates onto 
the disk. This is an important prerequisite for studying 
gravitational instability and frag mentation in young pro- 
toste l lar disks at all epochs (e.g. .IVorobvov fc Basull2006l 
[20101 iMachida et alJl2010t iKratter et al.ll2010D . 

We introduce a sink cell at r sc = 6.0 AU and impose 
a free inflow inner boundary condition. In the early pre- 
stellar phase of evolution, we monitor the mass accre- 
tion rate through the sink cell and introduce a central 
point-mass object (representing the forming star) when 
the mass accretion rate reaches a peak value (see Fig. |H] 
for details). In the subsequent evolution, approximately 
95% of the accreted material lands directly onto the star, 
while the rest is maintained in the sink cell in order to 
keep its density equal to the mean density of the gas in 
the innermost 1-2 AU of the actual numerical grid. The 
sink cell is otherwise dynamically inactive. It contributes 
only to the total gravitational potential, so that a smooth 
transition in column density is maintained from the nu- 
merical grid, into the sink cell, and to the protostcllar 
surface. 

We solve the mass and momentum transport equations, 
written in the thin-disk approximation as: 

<9£ 

_ + V p -(E Wp ) = ! (1) 



- (Zv p ) + [V • (E« p ® « p )] p = -V P V + £g p , (2) 

where the subscript p refers to the planar components 
(r, <p) in polar coordinates, £ is the mass surface den- 
sity, V = j Z z Pdz is the vertically integrated form of 
the gas pressure P, Z is the radially and azimuthally 
varying vertical scale height determined in each com- 
putational cell using an assumption of local hydrostatic 
equilibrium (IVorobvov fc Basul I2009D . v p = v r r + v^tj) 
is the velocity in the disk plane, g p = g r r + g^tfi is 
the gravitational acceleration in the disk plane, and 
V p = rd/dr+<f>r~ 1 d/d(f> is the gradient along the planar 
coordinates of the disk. The planar components of the 
di vergence of the symmetri c dyadic T,v p <S> v p are found 
in IVorobvov fc Basul (|2010f ) . The thin-disk approxima- 
tion is an excellent means of studying the evolution over 
many orbital periods and across a wide parameter space. 

The gravitational acceleration g p includes contribu- 
tions from the central point- mass star (once formed), 
from material in the sink cell (r < r sc ), and from the 
self-gravity of the circumstellar disk and envelope. The 
gravitational potential of the circumstellar disk and en- 
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Fig. 1. — Temperature evolutio n of zero-meta l licity gas derived 
from the one-zone calculations of Omukai et al. (2005, solid line), 
and the piecewise polytropic fit employed within the simulations 
discussed herein (dashed line); the filled circles mark the critical 
turning points of the fit specified in Table [T] 

vclope is found by solving the Poisson integral 

/r ut 
r'dr' 

x r - mm (3) 

■>° Jr' 2 + r 2 - 2rr' cos(0' - cj>) 

where r out is the radial position of the outer compu- 
tational boundary, or equivalcntly, the initial radius of 
a pre-stellar core. This integral is calculated using a 
fast Fourier transform technique which applies the two- 
dimensional Fourier convolution theorem for polar co- 
ordinates with a logarith mically-spaced radial grid (see 
iBinnev fc Tremaindll98"7i Sect. 2.8). 

Equations (p} and (|2|) are closed with a barotropic 
equation o f state, based o n the ID core collapse sim- 
ulations of I Omukai et al.l ([2005) that included the de- 
tailed chemical and thermal processes of the collapsing 
gas. Figure [1] shows the gas t empe rature versus density 
relation from lOmukai et all (|2005l ) for zero metallicity 
(solid line) and our piecewise fit (dashed line), with the 
transition points denoted by filled circles. This piecewise 
polytropic form can be expressed as follows 

fc-i 

Pk = c s v 7fc n Pop i+i > f ° r p-^-i <p< pc,k , (4) 

1=1 

where c s = y/RJT/Jx is the initial isothermal sound speed, 
T is the initial gas temperature, 1Z is the universal gas 
constant, and fi = 2.27 is the mean molecular weight of 
the primordial gas 4 . The value of the index fc distin- 
guishes the six individual components of the piecewise 

4 True for gas with volume density higher than 10 10 cm" 3 . As a 
rule, our model disks are characterized by volume densities higher 
than this value, except perhaps in the very outer parts at r > 
300 - 400 AU. 
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TABLE 1 

Parameters of the Barotropic Relation 



TABLE 2 
Model Parameters 



k 


Ti 


Pc,\ 
(g cm -3 ) 


(cm -3 ) 


1 


1.00 


8.20xl0~ 19 


2.16x10 s 


2 


1.14 


4.19X10- 16 


1.10x10 s 


3 


1.08 


4.50xl0- 1() 


1.18xl0 14 


1 


1.01 


1.32xl0~ 7 


3.47xl0 ie 


5 


1.12 


1.23xl0~ 3 


3.24X10 20 


6 


1.42 







form. We note that when k = 1 the product term is 
unity, and the pressure reduces to P 1 = c^p 11 . As our 
simulations evolve the effective surface mass density of 
the gas, the corresponding form of the barotropic rela- 
tion used in the code is 



fc-i 



r k = ci e 



2y7fc 



II ^ 



7«+i 



for S c < E < E c 



(5) 



(=i 



where the transition surface and volume mass density 
are related to one another through the instantaneous lo- 
cal scale height Z at each point in the disk via p Ct i = 
H Ci i/2Z. The transition points k, the associated mass 
and number volume densities at which these transitions 
occur, p c i and n c i, and the value of the various associ- 
ated polytrope indices are given in Table [T] Note that 
£ c ,o and £ c ,6 in equation (J5j arc formally equal to zero 
and infinity, respectively. 

The initial gas surface density E and angular velocity 
fi profiles for the primordial cores are similar to those 
that have been considered in the context of p resent-day 
star formation (jVorobvov k Basull2006l I201C 



E = 



roE 



^Jr 2 +■ 



o ('— ^ 




(6) 



(7) 



The radial pr ofile of E is an integr ated form of a Bonnor- 
Ebert sphere (|Dapp k B asu 2009() . while that of is the 
expected differential rotation profile to accompany ([6]) for 
a core contr acting from near-uniform initial conditions 
(jBasul fl997| ). The parameters Qq, Eo, and ro, are the 
central angular velocity, central gas surface density, and 
the radius of a central near-constant-density plateau, re- 
spectively. The latter is proportional to the Jeans length 
and is defined as 

VJc 2 

The parameter A determines the level of gravitational 
domination in the initial state. It is set to 1.5 for all 
models considered herein, and the initial gas temperature 
is set to 250 K (unless otherwise stated). We note that 
we have previously shown that the qualitative features of 
the protostellar accretion and disk evolution are weakly 
sensitive to th e specific profiles adopted for the initia l 
prcstcllar core (|Vorobvov k Basull2009t IVorobvovll26"T2t ). 
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So 
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n> 


M c 


P 


T 




g cm~ 2 


(km s _1 pc -1 ) 


(PC) 


(Mq) 


(io- 3 ) 


(K) 


ref. 


0.26 


0.75 


0.4 


176 


2.76 


250 


1 


0.26 


0.375 


0.4 


176 


0.69 


250 


2 


0.5 


0.75 


0.29 


179 


1.97 


350 



To set up a model prestellar core using equations ([6]) 
and |(7J), the values of Eo and Ho need to be specified. 
The former is found from equation by assuming a 
constant ratio of the core radius r out to the radius of the 
central plateau ro: r out /ro = 6. All model cores therefore 
possess a similarly truncated form. The parameter ro is 
chosen so as to form cores with mass on the order of 
200 Mq. This value is typical for primordial collapsing 
starless cores, as found in th e num erical hydrodynamics 
simulations of lYoshida et al.l ([2006D . 

The central angular velocity il of our model cores 
are chara cterized by the dimcnsionlcss parameter i] = 
VLlrl/c 2 s (iBasul Il997h . which is related to the ratio 
/3 = E ro % 1 1 Sgrav | of rotational to gravitational energy 
by /3 ps O.977. The value of the central angular velocity Qq 
is then calculated by choosing a value of (3 appropriate 
for primordial cores with s pin parameter a = \/~j3~ 0.05 
(jO'Shea k Norman 1 12007| ). Table El shows a list of pa- 
rameters for the models presented in this study. 

All models are run on a polar coordinate (r, cf>) grid 
with 512 x 512 spatial zones. The inner and outer bound- 
ary conditions are set to allow for free outflow from the 
computational domain. The radial points are logarith- 
mically spaced, allowing for better numerical resolution 
of the inner grid, where the disk forms and evolves. In 
the reference model, the innermost cell outside the cen- 
tral sink has a radius of 0.11 AU and the radial and 
azimuthal resolution are about 1.9 AU at a radius of 100 
AU. This resolution is sufficient to fulfil the Truelove cri- 
terion, which states that the local Jeans length must be 
resolv ed by at least four numerical cells (| Truelove et al.1 
[1999] ). 

Indeed, the Jeans length of a thin self-gravitating disk 
can be written as 

For a surface density of E ps 500—1000 g cm~ 2 and tem- 
perature TrjI.O — 1.5xl0 3 K, typical for our disks at 
r« 100 AU (see Figs. Q] and [7]), the corresponding Jeans 
length varies between i?j«35 — 90 AU and is resolved by 
roughly 17-45 grid zones in each direction (r, 4>). 

Finally, we want to emphasize that we do not resort to 
the use of sink particles of any sort except for the central 
sink cell. The fragments studied in this paper are fully 
self-gravitating protostellar embryos supported against 
gravity by pressure and rotation. 

3. RESULTS 

In this section we consider the time evolution of 
our reference model (see Table EJ), starting from the 
prestellar phase and ending once the mass of the cen- 
tral star reaches 45 M Q . Beyond this mass, the ef- 
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Fig. 2. — Comparison of the core properties in the reference 
model with those found using numerical hydrodynamics simula- 
tions of collapsing primordial mini-halos. The solid lines present 
the radial profiles of the enclosed mass (top-left panel) and gas 
volume density (bottom-left panel), and also the specific angu- 
lar momentum and radial velocity vs. enclosed mass (top-right 
and bottom-right panels, respectively). The plus sign s and cr osses 
are the correspondin g data taken from IClark et al.l d2011aT l and 
lYoshida et alj j2006D . respectively. The dashed line shows the 
r~ 2,2 profile, typical for the volume density in collapsing mini- 
halos. 
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Fig. 3. — Gas surface density map in the reference model showing 
the disk when its age is just 800 yr. The scale bar is in log g cm~ 2 . 
Multiple fragments are already evident in the disk at this early 
time. 

feet of stellar irradiati on may strongly affect our results 
(|Hosokawa et al.ll2011l) . 

3.1. Cloud core at the onset of the formation of the 
central protostar 

We start by comparing the properties of our collapsing 
core in the reference model just prior to the formation of 
the central star to those derived using three-dimensional 
numerical hydrodynamics simulations of collapsing pri- 



mordial mini-halos. The solid lines in Figure [2] show the 
enclosed mass M enc vs. radial distance r (top-left panel), 
the gas volume density n g vs. r (bottom- left panel), the 
specific angular momentum L vs. M enc (top-right panel), 
and the radial velocity i> r vs. M enc (bottom-right panel). 
The volume density was retrieved from the gas surface 
density S and the vertical scale height Z using the re- 
lation n g = S/ '(2 Z '^mh). We compar e our model pro- 
files wjtdVthos^rf OaAleF^i] ()2011al ) (plus signs) and 
and lYoshida et all (|2006l ) (crosses). The form of the core 
in the reference model is similar to those derived using 
3D simulations. For instance, the density profile follows 
closely the usual r -22 form (shown by the dashed line to 
guide the eye) . A steeper falloff in the core outer regions 
is cau sed by a finite core boundary ([Vorobvov fc Basul 
l2005bl ). However, notable differ ences are seen in the ac- 
tual values of n g , M enc , v T , and L. For instance, the gas 
density in o ur model is about a factor of several lower 
than that of lYoshida et al.1 ()2006l) and about an o rder of 
magnitude lower than that of IClark et al.1 (|2011a[ ). This 
is primarily caused by the fact that these authors consid- 
ered more massive cores than we did in our study. This 
is evident from the top-left panel showing the enclosed 
mass as a function of radius. The specific angular mo- 
mentu m in our model is a f actor of several smaller than 
that of[ci ark ct al.[ (j2011af ) but is approaching that of 
lYoshida et al.l ( 20061 ) in the inner region. We could have 
reconciled the numerically derived profiles by considering 
cores with a higher initial positive density perturbation A 
and/or higher angular momentum. This, however, would 
not alter our main results because such cores would pro- 
duce more massive and extended disks, whose properties 
would favor gravitational frag mentation even more (e.g. 
iVorobvov fc Ba^l2lMl2lna) . 

3.2. Formation and evolution of a primordial disk 

The disk in the reference model starts to form at about 
t = 1.0 kyr after the formation of the central protostar. 
This delay between the formation of the protostar and 
the disk could have been even shorter had the sink cell 
radius been smaller than 6 AU. The first episodes of frag- 
mentation occurred about 300-400 yr after the formation 
of the disk. Figure [3] shows the gas surface density (in 
log g cm -2 ) in the inner 200 x 200 box when the disk 
age is only 800 yr. The red circle in the coordinate cen- 
ter represents the sink cell. Several fragments have al- 
ready formed by this time. Thi s fragmentat i on time scale 
is longer than th at found by IClark et al.l ()2011al ) and 
iGreif et "ail (|2011l) . who used sink cells with a smaller 
radius of 1.5 AU, but is comparable t o the fragmenta- 
tion timescale found by ISmith et al.l (|2011| ) , who used 
larger sink cells (20 AU). We note that our core is ini- 
tially of lower density and angular momentum than in 
those studies, which may act to increase the disk frag- 
mentation timescale in our simulations. 

Figure [4] presents a series of images of the gas surface 
density (logarithmic in g cm~ 2 ) for the inner thousand 
AU 5 . The time elapsed since the formation of the central 
star is indicated in each panel. We note that the whole 
computational domain extends to roughly 80,000 AU (ex- 
actly, 0.4 pc) but that we focus on the innermost regions 

5 The top row has a spatial scale of 1000x1000 AU, while the 
bottom row has a scale of 1400x1400 AU. 
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Fig. 4. — Gas surface density maps of the circumstcllar disk in the reference model. The time elapsed since the formation of the central 
star is indicated in each panel. The scale bar is logarithmic in g cm -2 . The yellow dashed lines indicate radial cuts passing through several 
fragments that are used to calculate the radial gas surface density profiles in Fig. [jj The red dashed curve illustrates an azimuthal cut used 
to calculate the gravitational torques acting on a fragment in Fig. [JJ and the arrow points to fragment F4 in Fig. \S\ 

where the disk forms around the central protostar. The 
number of fragments present in the disk varies with time, 
indicating that they are not long lived, either migrating 
onto the star or being dispersed. In both cases, gravita- ~ 
tional torques from spiral arms and/or other fragments < 
are responsible. However, new fragmentation episodes o 
take place because the disk is constantly supplied with 3 
material from the parent core. Most of the fragments ^ 
form in the intermediate and outer disk regions, in -%-2oo\ 
agreement with numerical m odels of disk evolution in <z 
the present-day un iverse fe.g.. lStamatellos &: Whitworthl 
l200cl [Clark e 2009). 



40 kyr 



The recurrent character of disk fragmentation is evi- 
dent in Figure |5] in which we plot the number of frag- 
ments N{ present in the disk with time. Nf increases 
steadily to 10 by t = 25 kyr, and then suddenly decreases 
by t = 30 Myr to Nf = 5, suggesting that half of the frag- 
ments have migrated or been dispersed in just 5 kyr. This 
is followed by a period of intense fragmentation followed 
by another rapid depletion during the subsequent 10 kyr. 
During the most vigourous episode of disk fragmentation 
the number of fragments in the disk peaks at Nf = 14 
at t = 50 kyr. This fragmentation burst is followed by 
a deep minimum when the number of fragments drops 
to just Nf = 2. It is important to note that Nf never 
drops to zero, suggesting that some of the fragments may 
not migrate inward but instead stay at quasi-stable wide 
orbits or even migrate outwar d as also found i n rec ent 
studies bv lClark et all (|2011a|) and lGreif et"all (|2011l) . 



-400 -200 200 400 -400 -200 200 400 

Radial distance (AU) 

Fig. 5. — Gas surface density maps for the reference model shown 
at t = 40 kyr (left-hand panel) and t = 45 kyr (right-hand panel) 
after the formation of the central star. The contour lines delineate 
regions in which the Toomrc Q-parameter is lower than unity. The 
scale bar is in log g cm -2 . 

The consecutive increases and declines in the number 
of fragments suggests that the disk approaches a limit cy- 
cle during which periods of vigourous disk fragmentation 
are followed by periods of relatively weak (to no) frag- 
mentation. This limit cycle behaviour is possible if the 
typical migration time scale of the fragments is shorter 
than the characteristic time required for disk loading of 
material inf ailing from the collapsing parent core. The 
disk therefore loses mass (via inward migration of the 
fragments onto the star) faster than it can be replen- 
ished via accretion from the parent core. Following sig- 
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Fig. 6. — Number of fragments in the disk at a given time, as 
seen in the reference model. 

nificant mass loss, the disk stays relatively dormant until 
sufficient mass is accumulated to trigger another burst of 
fragmentation. 

For a mass infall rate onto the disk of Ma w 
10~ 3 M Q yr -1 and disk mass Afd 10 Mq (see Fig. [5]), 
the characteristic timescale for disk mass loading is 

ti = ^& « 10 kyr. (10) 

M d 

The migration timescale t m estimated for fragments in 
the reference model lies in the range of 10 3 — 7 x 10 4 yr, 
with most fragments having migration times of just a few 
thousand years (see Table [3] and equation [12]). Hence, 
the mean migration time is indeed shorter than the char- 
acteristic mass loading time. It is worth noting that t\ 
is in approximate agreement with the temporal variabil- 
ity of 5 — 15 kyr over which peaks in the distribution of 
fragment numbers occur in Figure [6] 

To illustrate the effect of gravitational fragmentation, 
regions in which the Toomre Q-parameter fall below 
unity are plotted in Figure [5] (traced in black) . The Q- 
parameter is defined as Cgfl/irGT,. Evidently, all frag- 
ments are characterized by Q < 1. In addition, some 
parts of the spiral arms are also characterized by Q < 1 
but do not yet show signs of fragmentation. Indeed, the 
number of fragments at t = 40 kyr (right-hand panel) 
is greater than at t = 45 kyr (left-hand panel), indicat- 
ing that new episodes of disk fragmentation took place 
during the intervening period. We note that the Q- 
paramcter in the innermost 50 AU is greater than unity, 
which explains the stability of this region to gravitational 
fragmentation. 

The left-hand column of Figure [7] presents the gas sur- 
face density distributions calculated along the radial cuts 
shown schematically by yellow dashed lines in the upper- 
left panel of Figure [1] (t = 30 kyr) . Each cut has an an- 
gular width of 15°, and is centered on the peak density of 
the fragment and and its surrounding minidisk. The gas 
surface density is averaged over the grid zones that are 
overlayed by the cut and have the same radial distance. 
The top, middle and bottom rows correspond to cuts 1, 
2, and 3, respectively. The positions of the fragments 
are indicated by arrows. The characteristic r slope is 
shown as a dotted line for convenience. Evidently, the ra- 



dial profile of the gas surface density scales as £ oc r , 
which is typical for self-gravitating dis ks in presen t-day 
star formation (e.g.. IVorobvov fc Basull2007l l2009t ). The 
typical gas volume density averaged over the azimuthal 
angle is a few x 10 14 cm" 3 at r = 10 AU and it drops 
to about a few xlO 11 cm" 3 beyond 100 AU. These val- 
ues are in reasonable agreement with three-dimensional 
numerical hydrodynamics sim ulations of disk for mation 
around primordial protostars ([Clark et al.ll2011a| ). 

The right-hand side column in Figure [7J presents the 
gas surface density distributions (solid lines) and nor- 
malized gravitational torques (dashed lines) calculated 
along the azimuthal cuts passing through each of the 
fragments. One such cut is shown by the red dashed 
curve in the upper-right panel of Figurc[4]for illustration. 
The width of the cuts is 10 grid zones, which, depend- 
ing on the radial position, translates into a radial width 
of 15-55 AU. The gravitational torque is calculated as 
r = —m(r, <j>) d$>/d<fi, where m(r, <fi) is the gas mass in a 
cell with polar coordinates (r, 0). The gas surface density 
and integrated torque are averaged over the grid zones 
that are overlayed by the cut and have the same radial 
distance. 

Figure [7J shows that the gas surface density in frag- 
ments is higher by about two orders of magnitude as 
compared to local disk values. The behavior of the nor- 
malized gravitational torque is interesting in that it ex- 
hibits a steep rise to a maximum positive value on the 
leading part of the fragment (greater azimuthal angles) 
and a deep drop to a minimum negative value on the 
trailing part (smaller azimuthal angles). The absolute 
value of r is several orders of magnitude lower every- 
where else. This specific form of r in profile is caused 
by tidal forces acting on the fragment from the rest of 
the disk and, in particular, from the spiral arms within 
which the fragments are usually nested. The trailing part 
of the arm (with respect to the fragment) exerts a neg- 
ative gravitational torque and pulls the fragment back, 
while the leading part of the arm exerts a positive torque 
and pulls the fragment forward in the direction of rota- 
tion. The resulting tidal force tends to shear apart the 
fragment, but the fragment's own self-gravity prevents it 
from dispersing. 

To calculate the properties of the fragments, we de- 
signed an algorithm that located the fragments in the 
disk and calculated their mass Mf, maximum surface 
density S max , physical and Hill radii, Rf and i?u, mass 
within the Hill radius Mh, radial distance of the frag- 
ment rf, and the integrated gravitational torque acting 
on the fragments rf. In particular, the Hill radius is 

where M* and Mf are the masses of the star and frag- 
ment, respectively. The radius of the fragment Rf is 
calculated from the known area occupied by the frag- 
ment, assuming a circular shape. Additional details of 
the tracking algorithm arc described in the Appendix. 

Figure [U presents zoomed-in gas surface density maps 
of four typical fragments. Fragments F1-F3 are named 
for the radial cuts denoted by the yellow dashed lines 
passing through the corresponding fragments in FigurelH 
Fragment F4 is marked in Figure 2] by the arrow. Frag- 
ment positions, as determined by our fragment-tracking 
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Fig. 7. — Left column: Radial gas surface density profiles cal- 
culated along the radial cuts passing through several individual 
fragments shown in Fig. 2] by yellow dashed lines. The number of 
the corresponding cut is indicated in each panel. Arrows point to 
the radial position of the fragments. The dotted lines show an r~ 1,5 
radial profile for comparison. Right column: The azimuthal pro- 
files of the gas surface density (solid lines) and normalized gravi- 
tational torque (dashed lines) calculated along the azimuthal cuts 
passing through the same fragments. One of these cuts is shown 
schematically in Fig. [4] by the red dashed curve. 

algorithm, are outlined in black, and their properties 
listed in Tabled 

All of the individual fragment radii Rf are smaller than 
their respective Hill radii Rn (see Fig. [8j, indicating that 
the fragments are gravitationally bound objects. Frag- 
ment masses Mf lie in the range 0.02 — 0.38 M Q . As the 
mass located within the Hill radius Mh is always greater 
than Mf , the excess material can be interpreted as a mini- 
disk. The presence of minidisks is most evident around 
fragments F1-F3 in Figure [51 Fragment F4 is interesting 
in that its radius is nearly equal to the Hill radius, which 
is a manifestation of tidal stripping; it is also the least 
massive of the four fragments. In general, the closer the 
fragment is located to the central protostar, the smaller 
its radius and mass. 

The orbital dynamics of the fragments depends on 
the sign of the integrated gravitational torque acting 
on them. As indicated from the data in Table G2 all 
of the fragments in Figure [5] migrate inward. We note 
that occasional outward migration of the fragments as 
a result of N-body-like gravitational scattering was also 
found in our models. The migration time scale can be 
estimated as follows. A (small) change in the angular 
momentum of a fragment (mass Mf) on a Kcplcrian or- 
bit, due to a (small) change in the orbital distance of the 
fragment dri can be written as dL = ^MfVfdrf, where 
tif = (GM^/rf) 1 / 2 . The migration velocity of the frag- 
ment is then 



j o dL 

ar { I — 



dt 



dt 



MfVf 



(12) 



torques t acting on the fragment, the characteristic mi- 
gration time is 

(13) 



2T 



Noticing that dL/dt = T, where the latter quantity 
is calculated as the sum of all individual gravitational 



The estimated migration time scale from each fragment 
is listed in Table |3j Fragment F4 and F2 have the short- 
est/longest migration time scales, of 10 3 yr and 7x 10 4 yr, 
respectively. We note that these values are instantaneous 
migration time scales, and that their migration patterns 
may change as more envelope material accretes on to 
the disk. Nevertheless, these migration time scales in- 
dicate that most of the fragments that are seen in the 
disk at t = 30 kyr will eventually have migrated onto 
the forming protostar by the end of our simulations, at 
£ = 55 kyr. The number of fragments in Figure [4] docs 
not steadily decline in time because subsequent genera- 
tions of disk fragmentation continue to be fuelled by mass 
loading from the infalling envelope (e.g., at t = 50 kyr). 

Ultimately, fragments that survive tidal stripping 
and/or shearing will pass through the inner boundary 
of the simulation, after which the evolution may branch 
into one of two possibilities. If the fragment survives tidal 
stripping and shearing, it may settle into a close orbit and 
form a low-mass companion close to the central star. Al- 
ternatively, the fragment may merge with or be tidally 
destroyed by the central star, releasing its gravitational 
energy in the form of a strong luminosity outburst. The 
latter possibility has been suggested by the recent nu - 
merical hydrodynamics simulation of iGreif et al.l (|2012l) . 
We consider this phenomenon in more detail below. 

3.3. Accretion and luminosity bursts 

Temporal evolution of the mass accretion rate in the 
reference model is shown in the top panel of Figure [HI 
from the disk onto the protostar (solid line), together 
with the smoothly varying accretion rate of material from 
the envelope onto the disk at 3000 AU (dashed line). The 
formation of the protostar is marked by the sharp peak 
in the mass accretion rate; we fix this point as t = 
and count time forward from here. The protostar then 
accretes material from the immediate surrounding enve- 
lope at the rate of a few times 10 -3 M yr -1 , growing 
to 6.0 M Q within approximately 2 kyr. At this point 
the accretion is temporarily halted by the formation of a 
quasi-Keplcrian disk. 

Accretion quickly resumes as the disk grows in mass 
and the influence of gravitational torques acting within 
the disk continues to redistribute matter and angular 
momentum. We calculate the mean accretion rate in 
1,000 year intervals following the formation of the quasi- 
Keplcrian disk, finding M to decline steadily from ap- 
proximately 10 -3 to 10 -4 Mq yr -1 over the course of 
the simulation. However, this smoothed background of 
accretion is punctuated by significant burst events during 
which clumps of typically ~0.03 M in size are accreted 
on time scales of less than 100 yr (e.g., fragment F4 in 
Figure [8j). A select few individual burst events involve 
the accretion of clumps of ~ 1.0 M . This results in 
effective mass accretion rates of ^-dO" 1 M yr" 1 , with 
the most extreme rates as high as 1.0 M Q yr -1 . Such 
large fluctuations in the accretion rate are the character- 
istic signature of gravitationally-induccd fragmentation 
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Fig. 8. — Zoom-in images of four typical fragments present in the disk at t = 30 kyr in Fig. [4] Fragments F1-F3 are named after the 
radial cuts that pass through the corresponding fragments, while fragment F4 is indicated by the arrow. The black contour lines delineate 
the fragments as determined by our fragment-tracking algorithm, while the yellow line outlines the Hill radius. The scale bar is logarithmic 
in g cm~ 2 . 
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within the disk, and subsequent accretion of the frag- 
ments onto the protostar. 

The cumulative effect of these two disparate modes of 
accretion are summarized by the mass growth curves in 
the middle panel of Figure [9] The protostar mass (solid 
line) grows rapidly over ~5 kyr during an initial phase of 
smooth accretion, continuing to increase steadily there- 
after. However, the burst mode of accretion is also evi- 
dent through abrupt increases in mass, of a few M Q at a 
time. These increases are typically followed by plateaus 
during which the mass of the protostar changes very lit- 
tle while the disk equilibrates back into a quasi-Kcplcrian 
state. Each burst event is also mirrored by a correspond- 
ing decrease in the total disk mass (dashed line) . Overall 
however, the disk mass continues to increase due to the 
accretion of material from the remnant parent cloud core. 

An object's luminosity during the early stages of pro- 
tostcllar evolution can be effectively characterized by the 
accretion luminosity alone. This is estimated by the dis- 
sipation of kinetic energy from infalling disk material 
landing on the protostellar surface: 



GAUM 
2R* ' 



(14) 



where M* is the mass of the protostar, M is the ac- 
cretion rate from the disk, and i?* is the protostellar ra- 
dius. To determine i?» in the absence of a detailed model 
for the stellar interi o r, we adopt the evoluti onary model 
of lOmukai fc Pallal ([200l . and following iSmith et al.l 
(|2011[) . employ simple power-law approximations for the 
protostar radius during each phase of its evolution. 

Following the initial c ollaps e of the cloud core, 
IStahler. Palla. fc Salpeterl (|1986f ) showed that the pro- 
tostellar radius grows according to a mixed power-law 
of the form R*ocM°- 27 M^ 1 ; for notational convenience 
we use M_ 3 to denote the ratio of the actual mass ac- 
cretion rate M to the fiducial value of 10 -3 Mq yr _1 
(|Omukai fc Pal la 2003). However, this growth is exac- 
erbated once the internal temperature of the protostar 
rises sufficiently to drive a wave of luminosity outward 
from the core. The result is a sudden and rapid expan- 
sion of the stellar surface. Once this wave reaches the 
surface of the star itself, the interior is able to relax, 
and the protostar begins Kelvin-Helmholtz contraction 
toward the main-sequence. The following power-law rela- 
tions therefore approximate the evolution of i?* through 
these distinct phases, as well as the transitions between 
them (jSmith et al.ll2011[) : 



R.„ 



26M°- 27 M°- 3 41 , 

AiMl 

AoMr 2 , 



Ah < Mi 

Mi < M* < M 2 

M 2 < M„ & i?* < R u 



(15) 



The constants Ai and A 2 are matching conditions that 
ensure the functional form of i?* remains smoothly vary- 
ing during the transitions. 

The mass parameter Mi marks the transition between 
the adiabatic phase of growth and the arrival of the lu- 
minosity wave at the protostellar surface; M 2 , the transi- 
tion between the luminosity wave driven expansion and 
subsequent Kelvin-Helmholtz contraction. Mi and M 2 
are fixed by the instantaneous mass accretion rate as the 
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Fig. 9. — Top: Mass accretion rates, from the disk onto the 
central protostar (solid), and from the remnant cloud core onto 
the disk at 3000 AU (dashed). Middle: Temporal evolution of 
the masses of the protostar and its disk. Punctuated increases in 
the protostar mass, associated with burst events, correspond to the 
momentary decreases in disk mass. Bottom: Accretion luminosity 
associated with the mass accretion rate onto the protostar. 

protostar transitions between phases, and are defined as 



Mi 
M 2 



5M° 3 27 , 
7M°_f. 



(16) 



Although this model was originally developed under 
the assumption of a constant mass accretion rate, evo- 
lution of the interior can be assumed to occur roughly 
adiabatically due to the long cooling time of the proto- 
stellar interior. As a result, significant rapid variability 
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in i?* is not expected while the Kelvin-Hclmholtz time 
scale is much longer than the accretion time scale. 

The bolometric accretion luminosity L acc for our refer- 
ence model is presented in the bottom panel of Figure [9l 
L acc reaches a value on the order of 10 4 Lq almost imme- 
diately at the time the protostar is formed, and continu- 
ing to vary slowly about this level over the course of the 
simulatio n, typical for proto stars of primordial composi- 
tion (e.g.. lSmith et al1l2012f ). After the disk is formed at 
t = 2 kyr, the accretion luminosity demonstrates a highly 
variable pattern of behaviour with high-intensity bursts 
(up to 10 7 Lq) superimposed on the lower baseline value 
(~ 10 4 Lq). This highly variable luminosity may have 
important consequences for the disk evolution, but more 
accurate numerical simulations taking into account the 
heating/cooling balance are needed to assess this effect. 

Spiral arms are formed when the destabilizing effects of 
self-gravity become comparable to, or dominate over, the 
stabilizing effects of pressure and shear. If the Toomre 
Q-parameter within the arm drops below unity, sections 
of the arm may collapse to form bound fragments. The 
accretion bursts occur when these fragments are driven 
onto the protostar by the gravitational torques of the spi- 
ral arms and/or other fragments. The temporal evolution 
of Q and the integrated torque T can be used to illustrate 
this phenomenon, with Q serving as an approximate sta- 
bility criterion, and T roughly expressing the efficiency 
of angular momentum and mass redi stribution by spi- 
ral in homogeneities within the disk (jVorobvov fc Basul 
120011 . As we are interested in the global properties of 
the disk (instead of their local variations), we calculate 
an approximate global value of Q = c s Q/ (nGT,), averag- 
ing c s , f2, and £ over all of the computational grid zones 
of the inner 500 AU region within which the disk is lo- 
calized. The quantity c s = (dP/cE) 1 / 2 is the effective 
sound speed. 

Figure [TTJ presents the mass accretion rate onto the 
protostar (in black) and the integrated gravitational 
torque normalized to a local maximum value (in red) 
as a function of the time elapsed since the formation 
of the protostar from our reference model. We focus 
on a short period of the evolution in order to make the 
interdependence between the bursts, integrated torque, 
and Q-parameter clearly visible. For the same reason, 
we also plot the mass accretion rate on a linear scale 
here. The red line shows that Q is smallest just before a 
burst, and its value rises sharply afterward, reflecting a 
transient decrease in the disk mass caused by fragments 
passing through the inner computational boundary. The 
integrated torque grows before each burst, reaching a 
maximum value at the moment of the burst itself due 
to strong gravitational interaction between spiral arms 
and fragments as the latter are torqued onto the pro- 
tostar. After the burst, the integrated torque returns 
to a marginal value. We note that the Q-parameter in 
Figure fTOl does not explicitly drop below unity — a typical 
value for gravitational fragmentation — because of the av- 
eraging over the entire disk; local values of Q (Figure [4} 
do drop below unity in the densest parts of the spiral 
arms and in the fragments. 

3.4. Parameter space study 

We next consider two additional models with the pur- 
pose of exploring the robustness of the accretion and lu- 




FlG. 10. — Top: Mass accretion rate onto the protostar (black 
solid line) and the global Toomre Q-parameter (red dashed line) 
versus time elapsed since the formation of the protostar. Bottom: 
Mass accretion rate onto the protostar (black solid line) and the 
integrated gravitational torque by absolute value (red dashed line) 
versus time. See text for an explanation of the interdependence 
between these quantities. 



minosity burst phenomena. We choose primordial cores 
with a lower initial angular momentum (Model 1) and an- 
other with a higher initial temperature (Model 2) than in 
our reference model. These choices are motivated by nu- 
merical studies of the burst phenomenon in present-day 
star formation that have shown both the number and in- 
tensity of such bursts to diminish with decreasing angu- 
lar m omentum and increasing temperature of the parent 
core ([Vorobvov fc Basul 120 10T ) . The specific parameters 
of each model are listed in Table [2j 

Figure QT] presents the mass accretion rates (top row), 
disk and stellar masses (middle row) , and accretion lumi- 
nosity (bottom row) for Model 1 (left column) and Model 
2 (right column). We terminate the simulations when 
the mass of the central objects reaches 45 M Q . Although 
the burst phenomenon is still occurrent, the strength of 
the accretion and luminosity bursts is somewhat lowered 
in comparison with our reference model. Though some 
of the luminosity bursts exceed 10 Lq in the reference 
model, in Models 1 and 2 most of the bursts stay between 
10 6 -10 7 L . 

The apparent decline in the strength of the bursts in 
Model 1 is caused by a notable decrease in the corre- 
sponding disk mass as compared to the reference model. 
Model 1 is characterized by a factor of four smaller ratio 
of the rotational to gravitational energy (3 than in the 
reference model, which leads to a delayed formation of 
the disk and overall lower disk mass. Low mass disks are 
less prone to gravitational fragmentation than their high 
mass counterparts, and the masses of their masses are 
correspondingly lower. The solid line in the middle-left 
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Fig. 11. — Mass accretion rates (top row), disk and stellar masses 
(middle row), and accretion luminosity (bottom row) versus time 
elapsed after the formation of the protostar for Model 1 (left col- 
umn) and Model 2 (right column). The solid lines in the top row 
show the mass accretion rate onto the protostar, while the dashed 
lines yield the mass infall rate at 3000 AU. The solid and dashed 
lines in the middle row are the masses of the protostar and the 
disk, respectively. 

panel in Figure [TT] indicates that the typical mass of the 
accreted fragments in Model 1 does not exceed 1.0 M Q , 
whereas in the reference model the fragments are of a 
few solar masses by the time they are accreted onto the 
protostar. This explains the somewhat weaker strength 
of the bursts in Model 1. 

Model 2 has a higher initial temperature (T = 350 K) 
and hence a higher disk temperature (because of the 
adopted polytropic equation of state) than the reference 
model (T = 250 K). For instance, the gas temperature 
at r = 100 AU in model 2 at t = 30 kyr is approximately 
1500 K, while in the reference model it is ps 1100 K. 
The contrast becomes stronger closer to the protostar. 
High temperature disks are also less prone to gravita- 
tional fragmentation because they have a higher Toomrc 
Q-parameter than their low temperature counterparts. 
This effect is somewhat offset by a higher mass infall 
rate onto the disk as indicated by the dashed lines in 
Figure 1111 As a result, the variations in accretion and 
luminosity in Model 2 are of smaller amplitude than in 
the reference model, but the burst phenomenon is still 
well pronounced. We also ran a model with a factor 
of 1.5 smaller initial core mass than that of the refer- 
ence model but similar /3, which yielded essentially sim- 
ilar results to those of Model 1. Obviously, models with 
increased /3 and/or M c that would yield an even more 
intense profile of bursts are not considered here. We 
conclude that the accretion and luminosity bursts are 
a robust phenomenon, which is expected to occur for a 
range of initial primordial core masses, angular momenta 
and temperatures. 

4. DISCUSSION 



Our parameter study successfully demonstrates the ro- 
bustness of the burst accretion phenomenon across a wide 
range of configurations in mass, rotation rates, and tem- 
peratures. We find the formation of a quasi-Kcplerian 
disk to be the outcome of the gravitational collapse of 
primordial prestellar cores. The disk expands quickly in 
radial extent to several hundred AU within just a few kyr 
of t he formation of the central protostar. In agreement 
with lClark et al.1 ()2011af ). as ha s been found in studies o f 
present day star formation (e.g.. lVorobvov fc Basul2007l) . 
we find that the disk exhibits a near-constant Toomre Q 
parameter ~ 1 . Figure [TTJ demonstrates that this is not 
only true initially, but holds over a substantial tract of 
time, ~80 kyr. 

Studies of mass accretion during primordial star for- 
mation have focused primarily on the smooth accretion 
of material onto the protostar, remini scent of Larson- 
Pens ton type an alytic solutions (e. g., lOmukai fc Pallal 
120031 ). Recently, ISmith et al.l (|2011[ ) studied the multi- 
plicity of fragmentation events in somewhat larger cloud 
cores. They found the accretion rates varied due to the 
motions of the protostellar fragments through regions of 
alternately high- and low-density gas within the natal en- 
vironment, with initial accretion rates of 10 -2 M Q yr _1 
that then declined to 10 -3 M© yr _1 ; in approximate 
agreement with analytic solutions. Indeed, we find simi- 
lar mass accretion rates onto the disk of ^10 -3 M© yr _1 . 
However, this source of infalling material is responsible 
for the mass loading of the disk and its subsequent frag- 
mentation (Section 3.1), which is the ultimate origin of 
the variability in the actual rate of accretion onto the 
protostar. 

The burst mode of accretion thus represents an en- 
tirely novel mechanism in contrast to either the clas- 
sical analytic scenario of monolithic accretion, or that 
due to the motion of individual protostellar fragments. 
The resupply of material infalling onto the disk at large 
radii establishes a pattern of recurrent gravitational in- 
stability driven fragmentation, as has been seen in mod- 
els of present-day disk evoluti on (e.g., IVo robvov fc Basul 
I2006D . In analogue studies by ISmith et al.l (|2012l ) it has 
been noted that the number of fragmentation events ac- 
tually increases steadily in time. We find a similar steady 
increase initially (Fig. [7]). However, SPH implementa- 
tions may be limited in their ability to resolve small-scale 
torques that may act to shear objects apart, and long- 
term temporal resolution is required to resolve the vari- 
ability and seeming periodicity in the rate of fragment 
formation observed in the simulations presented herein. 

The large-amplitude variations in the mass accretion 
rates are also responsible for giving rise to a corre- 
spondingly large-amplitude time varying accretion lu- 
minosity (Fig. [8]). Although the accretion luminosity 
associated with primordial star formati on is often re- 
garded as only a mild he ating source (jMcKee fc Tanl 
120081 : iHosokawa et al.ll201lT) . the nature of episodic accre- 
tion suggests that the accretion luminosity itself may ex- 
ceed prior estimates by a factor of between 10-100 times. 
This may raise the profile of the heating rate due to ac- 
cretion to be comparable to that from gas compression. 
However, our polytropic modeling of the gas thermody- 
namics is the primary caveat limiting our ability to an- 
alyze the complete effect of this heating mechanism on 
the disk's evolution. The robustness of our simulation re- 
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suits, with respect to higher disk temperatures however, 
suggest that this mechanism may be self-regulating, with 
increased heating due to bursts stymieing fragmentation, 
allowing the disk to cool, and in-turn increasing again the 
fragmentation likelihood. 

Although most fragments are torqued inward and onto 
the central protostar, contributing to the luminosity 
bursts, several alternative fates also exist. Fragments ca- 
pable of settling into stable orbits may evolve in conjunc- 
tion with their hosts to fo rm binary -pairs of firs t stars 
(e.g., iMachida et al] 12008(1 . In fact, ISuda et al.l (|2004f l 
has shown that the abundance patt erns of at least two hy- 
per metal-po or _QFe/H]<-5) stars (jChristlieb et al.ll200ll : 
iFrebel et al.ll200l can be explained by a unique history 
of mass transfer onto a first-generation low-mass binary 
star. 

Finally, we note that several of the fragments (F2 for 
example, in Figs. [4] and |8|) are seen being ejected to 
substantially larger orbits as a result of N-body-lik e in- 
teractions. Recen t studies by IClark et al.1 (|2011af l and 
iGreif et"aL1 (| 2 1 ll) noted similar dynamically chaotic in- 
teractions between fragments formed in those simula- 
tions. Although these fragments are most often sheared 
apart in the outer disk, we posit that a low-mass 
fragment — that is itself self-gravitating — may be ejected 
from the protostcllar embryo as seen in simulation s 
of present-day star formation (|Vorobvov fe Basv3l2012[ ). 
This raises the tantalizingly possibility of a heretofore 
unseen population of low-mass primordial stars that may 
have survived into the present-day. 

5. MODEL CAVEATS 

The thin-disk approximation. The applicability of the 
thin-disk approximation in our models is discussed in 
the Appendix. Here, we want to stress three additional 
points. First, the aspect ratio A of the disk scale height 
Z to radial distance r strongly depend on the disk-to- 
star mass ratio £ (see equation lAlj) and consequently 
on the initial conditions in the primordial cores. Mas- 
sive cores with high angular momentum are expected to 
form massive disks soon after the formation of the pro- 
tostar. These disks can be character ized by high £ ~ 1. 
and, consequently, by high A ~ 1.0 (|Clark et alj|2011ah . 
Our models have £ < 0.5, which allows us to use the 
thin-disk limit. Second, the disk structure is quite ir- 
regular and the thin-disk approximation may break lo- 
cally, even though it is fulfilled globally. And finally, full 
three-dimensional n umerical simulations of present-day 
star formation (e.g. IMachida et al.ll2011|l confirmed ro- 
bustness of the burst phenome non originally discovered 
using the thin-disk simulations (jVorobvov fc Basull2005aL 
1 2 6h . Three-dimensional simulations of primordial star 
formation already showed quic k inward migratio n of the 
fragments on short timescales (|Greif et al.ll2012|l and we 
await confirmation of the repetitive nature of the burst 
phenomenon on long timescales. 

Barotropic equation of state. In the present study, we 
approximated the thermal balance of the gas using a 
barotropic equation of state. During the collapse of a 
cloud core the major heating source comes from com- 
pression and the barotropic approximation is justified. 
However, once star formation is underway, radiative 
cooling and stellar irradiation may produce a range 
of temperatures that cannot be explained by a simple 



barotropic approximation (jClark et alj|20lil) . Numerical 
simulations of fragmenting barotropic disks seem to 
yield more fragments than tho se that take into account a 
detai led thermal balance (e.g. IStamatellos fc Whitworthl 
12009( 1. We note, however, that the burst phenomenon 
in the present-day star formation, ori ginally discovered 
using a barotropic equation of state (|Vorobvov fc Basul 
l2ooeh . was later shown to exist when more detailed 
thermal physics calcula tions were taken into account 
(|Vorobvov fc Basull2010ft . 

Stellar irradiation. A burgeoning protostar may start 
affecting its surroundings quite early in the evolution. 
For the mean accretion rate of 10~ 3 Mq yr _1 , a charac- 
teristic upper limit in our models, stellar UV irradiation 
becomes notable at a mass of about M * = 15 M & 
(iSmith et al.ll20Tl ISosokawa et al.ll20Tl . O ur numer- 
ical simulations extend to the time instant when the 
mass of the protostar reaches 45 Mq, meaning that the 
disk evolution may be affected by the lack of stellar 
irradiation. Numerical simulations taking into account a 
more detailed thermal balance are needed to access the 
effect of stellar irradiation on the burst phenomenon. 

6. CONCLUSIONS 

We have investigated the gravitationally-induced col- 
lapse of prcstcllar cores having pristine primordial 
gas composition, using nonaxisymmetric hydrodynam- 
ics simulations in the thin-disk limit. The gas thermal 
chem istry is modelled with a barotropic relation adapted 
from lOmukai fc Pa lla (2003j). We follow these simula- 
tions, in the absence of protostellar feedback, to the point 
at which UV ionizing irradiation from the central star 
would become important (while M < 45 M©). Our main 
conclusions are as follows: 

• Recurrent gravitational instability driven fragmen- 
tation and accretion of the fragments is an impor- 
tant mechanism through which protostars accumu- 
late mass in the early universe. This mechanism 
is mediated by smooth mass infall from the sur- 
rounding envelope at M« 10 -3 M Q yr _1 . As mass 
is loaded onto the disk, gravitational instability 
eventually induces fragments to form, and these 
are then driven onto the protostar — resulting in M 
of up to 10 _1 M Q yr~ : — by gravitational torques 
acting within the disk. 

• The burst mode of accretion is sensitive to varia- 
tions in the initial conditions (mass, rotation rate, 
and temperature) of the parent core, with the 
strength of bursts decreasing for decreasing (3 — the 
ratio of the magnitudes of rotational to gravita- 
tional energy — and/or increasing disk temperature 
T. 

• We have considered cloud core masses and rota- 
tion rates that are somewhat lower than those 
derived from the numerical hydrodynamics sim- 
ulations of collapsing primordial mini- halos (e.g. 
lYoshida et all 120061 : IClark et al.ll2011al) . Even in 
this case however, the burst phenomenon could not 
be entirely suppressed, and is expected to be more 
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Fig. Al. — Aspect ratio A of the disk vertical scale height to radius (Z/r) as a function of radius (r) in the reference model at t = 30 kyr 
after the formation of the central star. The solid line p rese nts the data calculated using the assumption of local vertical hydrostatic 
equilibrium and the dashed line is derived using equation I All . 



vigourous for other possible initial conditions. We 
conclude that the burst mode of accretion is likely 
a robust phenomenon in primordial star formation. 

Accretion luminosity produced by episodic mass 
accretion may contribute significantly to the heat- 
ing of material in the immediate vicinity of the pro- 
tostar, but numerical simulations taking into ac- 
count the heating/cooling balance are needed to 
assess this effect. 

Not all fragments eventually migrate inward and 
onto the protostar; some of the fragments are seen 
being ejected to the outer re gions of the disk, in 
line with previous stu dies bv IClark et all (|2011a| ) 
and iGreif et "ail (|2011[ ) . Given the implications for 



low-mass primordial star formation, the likelihood 
for survival of these fragments needs to be studied. 
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APPENDIX 
THIN-DISK APPROXIMATION 

The applicability of the thin-disk approximation hinges upon the geometry of the considered configurations. Elon- 
gated shapes are typical for protostcllar disks in nearby star-forming regions. Specifically, the thin-disk approximation 
is well justified as long as the aspect ratio A = Z/r of the disk vertical scale height Z to radial distance in the plane of 
the disk r does not consider ably exceed 0.1. This condition is usually fulfilled for disks with sizes up to 1000 AU (sec 
e.g.. IVorobvov fc Basull2010| ). but its applicability may be less obvious for the evolution phase immediately preceding 
disk formation, when material accretes directly from the infalling core onto the forming star. Nevertheless, the tem- 
poral behaviour and magnitud e of the accretion rate onto the star is very similar in the pre-disk stage for spherically 
symmetric and disk-like cores (fybrobvov fc BasvJl2005aL I2001I . ensuring that the stellar masses are calculated accu- 
rately in our simulations. Below we provide some analytical and numerical estimates justifying our use of the thin-disk 
approximation. 

In a Keplerian disk, the aspect ratio A = Z/r can be expressed as (|Vorobvov fc Basull2010l ) 



A < 



Qcr M d (r) 
CM* 



(Al) 



where Md(r) = J S(r, <j)) r dr d<f) is the disk mass contained within radius r, M* is the mass of the central star, Q cr 
is the critical Toomre parameter, and C is a constant, the actual value of which depends on the gas surface density 
distribution S in the disk. For a disk of constant surface density, C is equal to unity and for the S cx r -15 scaling 
typical for gravitationally unstable disks (see Fig. [7]), C = 4. Adopting a conservative value of C = 2 and setting Q cr 
to unity — characteristic for fragmenting disks — we obtain the radial profile of A shown in Figure IA1I by the dashed 
line. The ratio Md(r)/M* was calculated using the reference model at t = 30 kyr. Evidently, the aspect ratio derived 
using analytic considerations is smaller than 0.2 in the inner several thousand AU, which validates our use of the 
thin-disk approximation. We note however that this approximation may become only marginally valid at large r 
where Md(r)/iW* approaches its maximum value. 

This estimate is confirmed by the exact calculations of the disk scale height in our models. The azimuth-ally averaged 
radial distribution of the aspect ratio A = Z/r in the reference model at t = 30 kyr after the formation of the central 
star is shown by the solid line in Figure IAT1 The vertical scale height Z is calculated assuming local vertical hydrostatic 
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equilibrium in the disk using the method described in IVorobvov fc Basul (|2009t) . Figure IA1I reinforces our analytical 
estimate and demonstrates that the thin disk approximation is certainly obeyed within the actual disk, which, according 
to Figures 2] and [7J does not exceed 1000 AU in radius. Within this radial extent, the corresponding aspect ratio is 
below 0.25. Only at radial distances well in excess of 1000 AU might the thin-disk approximation be violated. 

Finally, we note that the thin-disk approximation assumes the absence of vertical motions, which turns the usual 
momentum equation for the z-componcnt of the gas velocity into the equation describing the local vertical hydrostatic 
equilibrium. This assumption is used to calculate the scale height Z and the aspect ratio A in equation (|A1|) . Therefore, 
a possible thick disk that is not in vertical hydrostatic equilibrium is not accessible through our modeling. 

FRAGMENT IDENTIFICATION ALGORITHM 

In the absence of sink particles, fragment identification on the computational mesh becomes a challenging task. 
Although from a numerical point of view the fragments are no different than from the rest of the disk, they can be 
identified based on the following set of physical criteria. First, we scan the disk and locate local maxima in the gas 
surface density that satisfy 

s> Mirar 5 < < B1 ) 

where Eioo is the gas surface density at a distance of 100 AU from the central star, and r is the radial distance in 
AU. A value of Eioo = 5000 g cm" 2 is chosen to represent typical densities of the fragments at 100 AU, based on the 
radial gas surface density profiles shown in Figure [7J Consequently, this criterion helps to filter out local maxima, e.g., 
spiral arms, while retaining those that represent the true fragments, which are usually characterized by a much higher 
density than the rest of the disk. 

After the radial and angular coordinates (r c , <j> c ) of the local maximum representing the center of a fragment have 
been identified on the computational mesh, we determine which of the neighbouring cells also belong to the fragment 
by imposing the following two conditions on the gas pressure P and gravitational potential $ 

dP 1 dP „ 

eP + ?W <0 ' (B2) 
a$ i 9$ . . 

d? + ?W >0 > (B3) 

where r' — r — r c and 4>' = <f> — 4> c . The first condition mandates that the fragment must be pressure supported, with a 
negative pressure gradient with respect to the center of the fragment. The second condition requires that the fragment 
be kept together by gravity, with the potential well being deepest at the center of the fragment. Although substantial 
support against gravity may be provided by rotation, we assume it to not invalidate the first criterion. 

In practice, we start from the center of the fragment and proceed in eight directions (along the coordinate directions 
and also at a median angle to them) until one of the above gradient criteria (or both) are violated. This procedure 
helps to outline an approximate shape of the fragment. We then check all the remaining cells that are encompassed 
by this shape and retain only those that meet both criteria. In addition, we filter out those cells with a gas surface 
density lower than that defined by equation (|B1[) with Si o = 5000 g cm~ 2 . 
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